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Abstract 

A method is proposed for determining the line tension, which is the main physical characteristic of 
a three-phase contact region, by Monte-Carlo (MC) simulations. The key idea of the proposed method 
is that if a three-phase equilibrium involves a three-phase contact region, the probability distribution of 
states of a system as a function of two order parameters depends not only on the surface tension, but 
also on the line tension. This probability distribution can be obtained as a normalized histogram by 
appropriate MC simulations, so one can use the combination of histogram analysis and finite-size scaling 
to study the properties of a three phase contact region. Every histogram and results extracted therefrom 
will depend on the size of the simulated system. Carrying out MC simulations for a series of system sizes 
and extrapolating the results, obtained from the corresponding series of histograms, to infinite size, one 
can determine the line tension of the three phase contact region and the interfacial tensions of all three 
interfaces (and hence the contact angles) in an infinite system. To illustrate the proposed method, it is 
applied to the three-dimensional ternary fluid mixture, in which molecular pairs of like species do not 
interact whereas those of unlike species interact as hard spheres. The simulated results are in agreement 
with expectations. 
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1 Introduction 

If two bulk phases are in equilibrium, the profiles of the order parameter (distinguishing between the phases) 
in the planar interface between them are non-uniform. Due to this non-uniformity, excess contributions to the 
extensive thermodynamic parameters of the system arise. These are usually referred to the dividing surface 
chosen within the interfacial region. An excess quantity per unit area of the dividing surface is called a surface 
adsorption. All surface adsorptions but one depend on the choice of the dividing surface. The special surface 
adsorption is that of the grand canonical potential, and is called the surface (or interfacial) tension; it does not 
depend on the choice of the dividing surface (in the case of a planar interface). 

The interfacial tension between coexisting phases plays a fundamental role in a theory of wetting phe- 
nomena and phase transitions. Various theoretical, numerical, and computational techniques exist for its 
calculation. Besides the mechanical definition of the surface tension, there are two statistical-mechanical 
recipes usually referred to as the virial route 1 and direct correlation function route (see, e.g., chapter 4 of 
ref.2). The latter is the more rigorous because it has no restrictions while the former is restricted to fluids 
with pair- wise additive intermolecular potentials. Although not easy in a general case, the equivalence of both 
methods can be proven (see, e.g., chapter 4 in ref.2 and references therein). These methods require the knowl- 
edge of the pair (and higher order) intermolecular potentials and the corresponding distribution functions in 
the interfacial region, which are not easy to obtain without drastic approximations. 2 

A more successful and traditional (although semi-numerical) way of calculating the surface tension is 
offered by density functional theory 2-4 which allows one to find the excess grand canonical potential asso- 
ciated with the interface between coexisting phases. This method has its own conceptual difficulties (see, 
e.g.,refs.3-5) related to (a) the choice of free-energy density for intermediate (between two bulk) values of 
the order parameter, particularly in the critical region where fluctuations are important, and (b) the long- 
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wavelength interface instability due to capillary waves. 

Finally, there is a variety of computer simulation techniques for calculating the surface tension. In the case 
of a system of continuous potentials the majority of simulations have used the mechanical definition 1 of the 
surface tension and the virial expressions for the normal and tangential components of the pressure tensor. 7-9 
The results obtained via these simulations show a large scatter. 10 Another simulation method involves the 
direct calculation of the free energy change in forming an interface. 11 ' 12 In simulation of lattice models a 
wider variety of methods has been used including thermodynamic integration using a biasing field to force 
the formation of the interface 13-15 and calculating the free energy differences for systems with periodic and 
anti-periodic boundary conditions. 16 ' 17 

All these methods implicitly or explicitly include conditions restricting the interfacial configuration space 
that can be sampled in simulations. Many of the methods also require establishing a well-defined interfacial 
region which can be quite difficult, particularly under conditions close to a critical point where the interface 
becomes very diffuse. An alternative approach to calculating the surface tension was proposed by Binder, 18 
whose method is based on the analysis of the probability distribution of states of a system as a function of 
the order parameter obtained from Monte-Carlo simulations as a normalized histogram. Under conditions 
where two phases can coexist in equilibrium, the histogram has two peaks each of which corresponds to a 
homophase state, i.e., to the whole system being in one of the two phases (Figure 1). Intermediate values of 
the order parameter correspond to heterophase states such that one part of the system is in one phase while 
the rest is in the other phase. The key idea of Binder's method consists of establishing a relation between the 
minimum of the histogram and the interfacial tension between two coexisting phases. A simulated histogram 
gives the interfacial tension of the finite-size (simulated) system, and the interfacial tension of the infinite 
system is obtained by extrapolating a series of results for systems of different sizes (but otherwise for the 
same thermodynamic conditions) to infinite size. Binder's method has been successfully applied to a wide 
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variety of models. 

If three phases are in equilibrium, there may exist a three-phase contact region in addition to three two- 
phase interfaces. The distortion of the order parameter(s) profile(s) in the three-phase contact region gives 
rise to excess contributions to the extensive thermodynamic characteristics of the system. These are usually 
referred to the three-phase contact line". Excess quantities per unit length of the contact line are called linear 
adsorptions. For three bulk phases in equilibrium, the interfaces between them can be treated as planar and 
the line of their intersection, that is, the contact line, as a straight line, so that its spatial location is determined 
by two coordinates in the plane perpendicular to it. The location of the contact line is arbitrary and among all 
linear adsorptions only one does not depend on its choice. This is the linear adsorption of the grand canonical 
potential also called the line tension. 2 All other adsorptions depend on the location of the contact line and 
hence are functions of its two coordinates. It should be noted that line tension may be positive or negative, 2 
unlike surface tension that is always positive. 

The line tension has been studied to a much lesser extent than the surface tension both experimentally 
and theoretically. 32 ~ 34 ' 2,35 ~ 40 There are no general statistical mechanical expressions for calculating the line 
tension that would be equivalent to the virial and direct correlation function routes in a theory of surface 
tension. For the case of two fluid phases in contact with a solid surface, a statistical mechanical theory 
for the line tension was proposed by Tarazona and Navascues; 41 using several assumptions about molecular 
interactions, they derived an expression for the line tension involving the one- and two-molecule distribution 
functions. All other microscopic models for the line tension were based on the local density approximation of 
density functional theory 2 ' 36-40 and are hence subject to conceptual problems similar to those for the surface 
tension. 3-5 As for computer simulation methods, they have often been applied to three-phase equilibria and 
wetting phenomena, 42-46 but we are not aware of direct molecular simulations of the line tension, which is 
probably partly due to the lack of a general statistical-mechanical theory thereof. It should be noted, however, 
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that studying a solid spherical particle at a liquid-vapor interface by molecular dynamics (MD) simulations 
and thus finding the contact angle and all interfacial tensions involved, Bresme and Quirke 43 evaluated the 
(particle-liquid- vapor) line tension with the help of a phenomenological expression for the particle free energy 
(found also by MD simulations) and modified Young's equation. 

In the present paper a method is proposed for calculating the line tension by MC simulations. The method 
is based on the statistical-mechanical analysis of the histogram of a system under conditions of three phase 
equilibrium involving a three-phase contact region. The paper is structured as follows. Section 2 outlines the 
original histogram analysis method for determining the interfacial tension of a two-phase interface pioneered 
by Binder. 18 Section 3 presents a general idea of how that method can be extended to three-phase equilibria in 
order to determine the line tension. Section 4 discusses a general set-up of MC simulation cell and boundary 
conditions and there is also developed a theoretical basis for the concrete realization of the key idea in a 
three-dimensional ternary fluid mixture, in which molecular pairs of like species do not interact while those of 
unlike species interact as hard spheres; 47-49 in this model a three-phase contact region may exist at densities 
above the "quadruple point" density. 47-49 The details of MC simulations and results for the surface and 
line tensions associated with the three-phase contact region are also presented in section 4. The results are 
discussed and conclusions are summarized in section 5. 

2 Binder's approach to determining the surface tension by Monte-Carlo sim- 
ulations 

Binder's method 18 for determining the interfacial tension of two coexisting phases is based on the order- 
parameter distribution formalism. Rather than set up an explicit interface in the MC-simulated system, this 
approach relies on spontaneous fluctuations that give rise to density inhomogeneities, which provide infor- 
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mation about the interfacial properties. The method has its own implementation difficulties, ' especially 
far away from the critical point, but they can all be alleviated by using special sampling techniques. 23 ' 31 ' 50 

Consider a finite "cube-shaped" system (to be simulated) of volume L d (d is the dimensionality) with 
periodic boundary conditions (PBC), which are necessary in order to avoid additional surface effects that 
can arise from the system boundaries. Under appropriate thermodynamic conditions (e.g., below the critical 
temperature if the temperature is a relevant thermodynamic parameter) the system will phase separate into 
phases a and (3 with the values m a and trip of the order parameter m. In a liquid-gas system the order 
parameter is the density; in a binary mixture of species A and B with a miscibility gap the order parameter is 
the mole fraction of one of the components in the coexisting A-rich and B-nch phases; in an Ising model, the 
order parameter is the spontaneous magnetization (two orientations thereof corresponding to the two phases). 
In all cases some conjugate field variable (the chemical potential, or the difference in chemical potentials of 
two components, or the magnetic field, respectively) has to be fixed at such a value that would correspond to 
coexisting phases in the thermodynamic limit. 

According to the theory of equilibrium fluctuations, 51 if the system is in a pure phase with the average 
order parameter m, the probability p{m) to find the system (with L — > oo) in a state with the order parameter 
m is given by the Gaussian distribution, 



where k is the Boltzmann constant, T is the temperature, and x is the partial derivative of the order parameter 
with respect to its conjugate field variable, with all other variables of state of the system being fixed. For the 
system at coexistence where both pure phases are equally likely to occur, the probability distribution of the 
order parameter can be represented as the sum of two displaced Gaussians, 



p(m) = L d / 2 (27tkT X y 1/2 exp[-{' 



m-m) 2 L d /(2kT X )}, 



(1) 



p{m) = iL d / 2 (27rfcr XQ )- 1 / 2 exp[-( 



m - m a ) 2 L d /(2kT Xa )] 
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+ h d / 2 (27rkTxpy 1/2 expl-(m-m(i) 2 L d /(2kTxp)}, (2) 

where Xa and x/3 are the susceptibilities in the pure phases a and (3. It is clear from eq.(2) that the probability 
of a homogeneous state with the order parameter m mm = (m a + mp)/2 decreases exponentially with the 
volume of L d of the system. On the other hand, the probability of heterogeneous fluctuations, such that two 
phases coexist in the system, decreases exponentially 18 ' 56 with the interface area 2L d ~ 1 , 

P« n ) oc p(m ) eM-^L^/kT], (3) 

where the superscript "ht" stands for "heterogeneous", the subscript "o" indicates either of the pure phases a 
or (3, and a is the interfacial tension in an infinite system. The size dependence of the pre-exponential factor 
p(m ) ~ \L d l 2 (2itkTxo)~ l I 2 is much weaker than the exponential one. Thus, for large enough L, the 
probability of "homophase" fluctuations with intermediate values m mm of the order parameter is negligible in 
comparison with the probability of "heterophase" fluctuations. 

As pointed out by Binder, 18 the size dependence of the pre-exponential factor in eq.(3) can be more 
complicated than that contained in p(m ), due to finite-size contributions to the interfacial free energy arising, 
e.g., from the "Goldstone modes" and capillary waves. As a result, one can expect that for large L the 
probability distribution p(m min ) (the superscript "ht" is omitted henceforth) has the general form 

p(m mm ) = aL x p{m ) ex^-toL^ 1 /KT\ (4) 
with unknown a and x, whence it follows that 



1 p(m ) a InL 1 



2L d-i ln kT + b L d-i + c L d-i > ( 5 ) 



(with unknown b and c) and 



—— = lim - , z ln — -. (6) 

kT L^oo 2L d - 1 p(m rain ) 
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This limit cannot be found analytically, but the probability distribution as a function of the order parameter, 
p{m ), can be obtained for a series of L's via appropriate MC simulations. Then, one can plot the LHS of 
eq.(5) vs lnL/L d_1 and the intercept of this plot with the ordinate axis will provide the ratio a/kT for the 
interfacial tension in the infinite system. 

3 Three-phase equilibrium: determination of the surface and line tensions 
by the histogram analysis method 

If the system is under such conditions that three phases can coexist in equilibrium, there are necessarily two 
order parameters such that at least one in the pair differs in different phases. For example, in a single compo- 
nent system at the triple point the gas, liquid, and solid phases can coexist, and a pair of order parameters can 
be the density and a structure parameter. In some two component systems ("solvent-solute") under appropri- 
ate conditions one can observe the coexistence of three fluid phases: solvent liquid, solvent vapor, and solute 
liquid. 30 In such systems the densities of the components serve as a pair of order parameters. For such sys- 
tems, the probability distribution of states is a function of two independent order parameters and determines 
a three-dimensional surface, hereafter referred to as a histogram (as in a single order parameter case). This 
three-dimensional surface has three peaks of equal height corresponding to three pure phases. 

Consider a three component fluid system of species 1, 2, and 3. Under appropriate conditions, such 
a system can phase-separate (demix) into three equilibrium phases differing from each other by a pair of 
independent order parameters, e.g., the mole fractions of any two components, say, m\ and ni2 (the third 
mole fraction is not an independent variable, m% = 1 — m\ — mi). Denote by p(m\,m<2) the probability 
distribution of states of this system with respect to the order parameters m\ and ni2. Figure 2a presents 
an example of the histogram (three dimensional surface determined by the function p(mi,m2)) generated 
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by semigrand canonical MC simulations in a cubic box with PBC for the three-dimensional ternary fluid 
mixture, in which molecular pairs of like species do not interact and those of unlike species interact as hard 
spheres (for details see section 4). Figure 2b shows the projection of the same surface onto the mole-fraction 
plane. The peaks of the histogram correspond to the pure phases, say, a, (3, and 7. The projections of the 
peaks are shown as the points 0\; (01 = a,/3, 7) in Fig.2b. Their coordinates m" 1 (i = 1,2) determine 
the composition (i.e. mole fractions) of the system in pure phases a, [3, and 7. If all the three phases are 
symmetric (as, e.g., for the model represented by the histogram in Fig.2), 

777" = = m], = = = rrig = mj = 777^, 

but this may not be the case in general. 

Suppose the system is in some pure phase o\ {o\ = a, /3, 7). According to the thermodynamic theory 
of equilibrium fluctuations, 51 the probability distribution of homogeneous fluctuations of order parameters 

mi,m2 is 



p( mi ,m 2 ) = v /det(re)— ^— exp 



2nkT/p 



V 2 
^11(777.1 - 771°) 



2kT/p 



+ 2Ki 2 (mi - 777°) (7772 - 777,2) + ^22 ("7-2 - ^2) ) 



(7) 



where V = L d is the volume of the system. In a semigrand canonical ensemble (SGCE), which we will use 
for our MC simulations, n is a 2 x 2 diagonal matrix with elements 



/ <9?77 7 ; 

= - 



(8) 



with A/jj = Pi — ^3 {i = 1,2) and pi the chemical potential of component i. The prime in eq.(8) 
indicates that one of the A/j's is kept constant in the partial derivative. The choice of SGCE is the matter 
of convenience and efficiency of MC simulations of phase separation in dense fluids 52-54 (see subsection 
4.2 for more details). Note that, in the thermodynamic limit, SGCE simulation at fixed V, T, N, Api,Ap 2 is 



equivalent to grand canonical simulation at fixed V, T, and such fj,i, [12, A*3 (fixed) that N, the grand-canonical 
average number of molecules in the system, is equal to N. 

As clear from eq.(7), the probability of a homogeneous state with one or both order parameters having 
their intermediate values, mi„ = (mf + mJ)/2 and ni2 mm = (mf + ml)/ 2, decreases exponentially with 
the system volume. Hence, for large enough volumes, "homophase" fluctuations with mi = m\ mm or/and 
m2 = ?ri2min are negligible compared with "heterophase" fluctuations. To show this, we will consider the 
system (to be simulated) that has the form of a cylinder of length L and radius R (Figure 3). Such a shape will 
be most convenient (but not necessary) to determine the line tension by computer simulations, but otherwise 
it does not affect measurable physical quantities because they are obtained by finite-size scaling, 18,19,31,53 " 55 
which involves their extrapolation to infinite cell sizes. 

3.1 Determination of the surface tension 

For "heterophase" fluctuations, intermediate points on the lines (a(3), ((3j), and (7a) in Figure 2b represent 
the two-phase states of the system with none of the third equilibrium phase present. For example, for a point 
on the line (a/3) phases a and (3 coexist while the amount of phase 7 is 0, although thermodynamically it 
can coexist with the a and (3. The middle points o a p = {mi mm , m2„), = {m,J,m2 mm }, and o ai = 
{mi„, ml} on these three lines correspond to the situations when one half of the system is in one phase and 
the other half in another phase. For example, the point o a p corresponds to such s situation that one half of 
the system is in phase a and the other half in phase [3, while the equilibrium phase 7 is not present at all. 
Let us use the generic symbol 02 to denote the points o a p, opj, o ia . The difference between the free energies 
of the system in states o\ and 02 is due only to an interfacial contribution F s from the interface between 
two coexisting phases in the state 02. According to the principles of equilibrium thermodynamics, p(o2) and 
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p(o±), the probability densities of two-phase states 02 and pure phase states o±, are related as 



p(o2)=p(o 1 )e- F °/ kT . 



(9) 



(phase 0\ is one of two phases in the state 02) 

There is no exact analog of PBC for simulations in a cylindrical geometry, but there are other boundary 
conditions attenuating the effect of cell boundaries. This issue will be discussed in section 4, but for now it 
need only be noted that one can impose PBC along the cylinder axis (hereafter called the z-axis) and assume 
that the effects of the lateral cell boundaries can be neglected. Because of PBC along the z-axis, there may be 
only three principally different "half-and-half spatial distributions of two equilibrium phases in the cylinder 
(Figure 3). Comparing the interfacial free energy contributions F s in these three cases, one can show that if 
the aspect ratio uj = L/R of the cylinder satisfies the inequalities 



then the phase partitioning shown in Fig. 3a is thermodynamically most favorable. Hereafter, this will be 
referred to as the axial "longitudinal" two-phase (AL2P) partitioning. The aspect ratio uj < 1 favors the 
"droplet-like" phase partitioning (Fig.3b), while uj > ir favors the "transversal" phase partitioning (Fig. 3c). 
For the droplet-like phase partitioning the interfacial tension is assumed to be the same as that of a planar 
interface. This should be a reasonable assumption for a large enough droplet. 

If the histogram is obtained by simulations in a cylinder subject to conditions (10) and PBC along the 
z-axis, one can represent eq.(9) in the form analogous to eq.(4), 



Here, a Q2 is the surface tension of a flat interface in an infinite system, and the pre-exponential factor aL x R y 
(with unknown a, x, y) is due to finite size effects (Goldstone modes, capillary waves, etc 18 ) that may con- 



1 < UJ < IT, 



(10) 



p(o 2 ) = aL x R y p{o 1 )e- { ? RL ' J °Jl kT . 



(11) 
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tribute to F s , where the leading term 2RLa Q2 is larger than finite size corrections by an order of magnitude. 
From eq.(l 1) it follows that 

I p(oi) a Q2 InL InR In a 

In — — - = — — — x — y — — . (12) 



2RL p{o 2 ) kT 2RL a 2RL 2RL 
Again, for an infinitely large system the RHS of this equation will equal the interfacial tension. The limit 
L — > oo, R — ► oo can be reached in such a way that u>(= L/R) remains constant, so that equation (12) can 
be rewritten in two forms, 

p(oi) a Q2 i JnL i 1 



2^L ln p(o 2 ) kT +b +C L 2 

°° 2 MR , 1 

= ir +d lP (13) 

with unknown 6, c, d, /. The probability distribution as a function of the order parameters, p(mi, m-i), can 
be obtained for a series of L's (or ii's) via appropriate MC simulations. Then, one can plot the LHS of 
eq.(13) vs InL/L 2 (or vs In R/R 2 ) and the ordinate intercept of this plot will provide the ratio a 02 /kT for 
the interfacial tension in an infinite system. 

3.2 Determination of the line tension 

"Heterophase" states, corresponding to points lying off the two-phase lines (a(3), (J3j), and (7a) in Figure 2b, 
represent three-phase states such that all three equilibrium phases are simultaneously present in the system. 
The amount of each phase present is determined by the coordinates of a point 03 corresponding to a given 
three-phase state. When the three equilibrium phases are symmetric, the central point e% of the composition 
triangle corresponds to a state in which all three phases are present in equal volumes. In a general case this 
three-phase volume equi-partition point is not necessarily in the middle of the composition triangle. 

Since the subject of the present work is the line tension, we restrict our interest to three-phase equilibria 
involving three-phase contact regions. The necessary and sufficient condition that there exists such a region 
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in a three-phase system is 

a ay < a af 3 + crg 7 , (14) 

(hereafter double Greek subscripts and superscripts indicate quantities for corresponding two-phase inter- 
faces). Condition (14) means that it is thermodynamically unfavorable for the system to have three phases 
arranged in three layers parallel to each other. When this condition is fulfilled, F s l, the difference between 
the free energies of the system in states o\ and o 3 may contain not only interfacial contributions from three 
two-phase interfaces, but also the "line tension" contribution from the three-phase contact region. Thus, p(o' 3 ) 
and p(oi), the probability densities of three-phase states o 3 and pure phase states o\ are related as 

p(o'3)=p( 0l ) e -^/ fcT , (15) 

where the excess free energy F s l has the form 

F sL = o afi A^ + afrAfr + a ai A^ + tL^\ (16) 

Here the ^4's are the surface areas of interfaces and L a ^ is the length of the three phase contact line (assumed 
to be straight) with which the line tension r is associated. 

Again consider a cylindrical simulation cell with the boundary conditions discussed above, and assume 
that there exists 1 < ojq < ir such that if L/R = ujq, then thermodynamically the most favorable three-phase 
partitioning is the "longitudinal" one shown in Figure 4a. When the three-phase contact line coincides with 
the axis of the cylinder (as shown in Fig.4a), this will be referred to as the axial longitudinal three-phase 
(AL3P) partitioning. Denoting by 03 the corresponding point in the plane m\,rri2, one can rewrite eq.(15) in 
the form 

p(o 3 ) = aL x R y p(o 1 )e-^ +CT ^ +C7 ^ LR+TL ^ kT , (17) 

where now cr's and r are the interfacial and line tensions, respectively, in an infinite systems, and the pre- 
exponential factor aL x R y (with unknown a, x, y, in general different from those in eq.(ll)) is again due to 
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finite size effects: because cr's and r in eqs.(16) and (17) are those of an infinitely large system, there may 
arise small (in comparison with the leading terms) finite-size corrections to F s l. 

Henceforth assuming that the simulations are carried out under such conditions that the AL2P and AL3P 
partitionings are thermodynamically most favorable, one can combine eqs.(ll) and (17) and obtain 

p(o 3 ) = aL x Ry\p{o 1 )\ 1 -Z\p{o 2 )fe- TL / kT , (18) 

where a, x, y are unknown (different from those in eqs.(ll) and (17)) and 

= cr af 3 + 0737 + o- a7 
2a Q2 

From eq.(18) it follows that 

For an infinitely large system the RHS of eq.(20) will equal the ratio of the line tension r to kT. If the limit 
L — ► oo, R — > oo is reached in such a way that uj = L/R remains constant, then eq.(20) can be rewritten as 

lln kit)]* = ^ InL 1 
L \p( 0l )]Z-ip(o 3 ) kT + L L 

with unknown b, c, d, f. Again, p(mi, m-i), the probability distribution as a function of the order parameters, 
can be obtained for a series of L's (or R's) via appropriate Monte-Carlo simulations at constant u(= L/R). 
Then, one can plot the LHS of eq.(19) vs lnL/L (or vs In R/R) and the ordinate intercept of this plot will 
provide the ratio r / kT for the line tension in the infinite system. 

It is worth emphasizing that the point o\ in eqs.(9)-(21) corresponds to any pure phase state (a, (3, or 
7) because (ideally) in a three-phase equilibrium ^(mf^mf) = pirn^.m^) = p(mj,m^). The point 02 
represents any two-phase equilibrium with the AL2P partitioning (Fig. 3a) involving o\ phase, and a 02 in 
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eqs.(ll)-(13),(19) is the corresponding interfacial tension. Finally, the point 03 in eqs.(17)-(21) represents 
the AL3P partitioning (Fig.4a). 



4 General set-up of the simulation cell and boundary conditions 

For simulations in a rectangular geometry one usually uses PBC in order to decrease the effect from the 
simulation cell boundaries. For simulations in a cylindrical geometry, there is no exact analog of PBS, but 
one can impose boundary conditions that do attenuate the effect of cell boundaries. Namely, one can impose 
the usual PBC along the z-axis (Figures 3 and 4). In the plane of polar coordinates 4>, r, perpendicular to the 
z-axis, impose "random-periodic" boundary conditions (RPBC); that is, if upon an elementary Monte-Carlo 
move the new coordinates z, (f>, r of a molecule are such that r > R (i.e., the molecule crosses through the 
lateral surface of the cylinder), the angle 4> is replaced by a new, random angle 4>' and the radius r is replaced 
by a new, periodic value r' = R— (r — R). After these changes are made, the acceptance criterion is checked 
and the move is either accepted or rejected. If the move is rejected, the molecule keeps its old coordinates 
(those prior to generating new coordinates z, cf>, r). In this work only hard-sphere type molecular interactions 
are considered, but in a general case (e.g., Lennard- Jones intermolecular potential) it would be necessary to 
deal with the problem of accounting for the interactions of "random-periodic" images of molecules, that are 
in the lateral surface layer of "cutoff thickness, with other molecules in this layer. 

The propensity of an equilibrium system to minimize its free energy and PBC along the z-axis prevent 
the formation of artificial (i.e., due to the finite size of the simulation cell) phase interfaces at the bases of 
the cylinder. However, imposing RPBC in the cf>, r plane gives rise to the probabilistic phase interfaces at the 
lateral surface of the cylinder. Actually, with RPBC imposed, a molecule crossing the lateral interface will 
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leave phase si and enter phase s 2 with the probability g slS2 (si, s 2 = a, (3, 7), 

g SlS2 = 9s 2S1 = ( s i> s 2 = a,/3,7), (22) 

(where si and S2 are not subscripts, they must be understood as the corresponding contact angles) with the 
normalization condition 9sxs 2 = 1 naturally satisfied. Thus, the whole lateral surface of area 2-itRL is an 
a[3 interface with the probability 2g a p, a (3j interface with the probability 2g@ 1 , and an aj interface with the 
probability 2g ai . Clearly, for any point o 2 (half-and-half two-phase coexistence) the analogous probability is 
2g 02 = 1/2. 

Thus, as an artifact of RPBC imposed at the lateral surface of the simulation cell (cylinder), it is necessary 
to modify equations (11)-(13), (17), and (19). Namely, eqs.(ll)-(13) must be replaced by 

p(o 2 ) = aL x i?^(oi)e-( 2+ ^) i?LCT °2/ fcT , (23) 

1 , p(oi) <Jo2 hi£ mi? In a 

m — — - = — — — x— — — — y— — — — — 7777" ) (24) 



(2 + ir)RL p(o 2 ) kT (2 + w)RL y (2 + n)RL (2 + tt)RV 
1 , p(oi) cr° 2 InL 1 



(2 + ir)RL p(o 2 ) kT L 2 L 2 

- + + (25) 

" kT + d R 2 + ; R 2 ' (jK>) 

respectively. Equations (17) and (19) must be replaced by 

p(o 3 ) = a^Ryp^y-^ 1 ^ 9 ^^^ 1 ^ 9 ^^^ (26) 
t _ (1 + 4:Trg a p)a af 3 + (1 + Ang^a^ + (1 + 47rg a7 )a a7 

(2 + 7T)(T 02 

Now, the ratio a , 2 / kT for the surface tension of an infinite system is determined as the ordinate intercept of 
the plot of the LHS side of eq.(25) vs InL/L 2 (or vs InR/R 2 ). The ratio r/kT for the line tension in an 
infinite system is found, as previously, by plotting the LHS of eq.(21) vs In L/L (or vs In R/R) and finding 
the ordinate intercept of this plot, but the parameter £ is now given by eq.(27) rather than eq.(19). 

16 



Note that even with perfect PBC imposed on the finite simulation cell, there may be cell boundary con- 
tributions to F s and F s l, but they can be expected to be much smaller than the above contributions from 
the physical phase interfaces at the cell boundaries. Moreover, in the case where the equilibrium phases are 
symmetric (such as Ising-like models, the model discussed below, etc), the cell boundary contributions to F s 
(and F s l) can be assumed to be zero 18 ' 55 because in such systems the free energy of the cell is affected by its 
boundaries in the same way, independently of which phases are present. 

4.1 Method specifics for the simulated model 

To illustrate the proposed method with a concrete example, MC simulations were carried out for the three- 
dimensional three-component fluid mixture, in which the molecules of different species interact as hard 
spheres, whereas molecules of the same species do not interact at all (i.e., are completely penetrable by 
one another). 47 ~ 49 This is a natural generalization of the binary model, earlier proposed by Widom and 
Rowlinson 56 and come to be known as the primitive version of the Widom-Rowlinson model. The inter- 
molecular potentials between two molecules of species i and j at a distance r from each other are defined 
as 

u u (r)=0 (i = 1,2,3), (28) 

u ij (r) = (r>r),i^j = 1,2,3), u tj {r) = oo (r < v , % + j = 1, 2, 3), (29) 

so that the molecules of each pure species constitute an ideal gas, while every molecule of species i sees 
every molecule of species j (j ^ i) as a hard sphere of diameter rf. The phase behavior of this system 
was analytically studied in both mean-field and higher order approximations. 47-49 The system is predicted 
to have a so-called "quadruple point" density p q : at this density within the composition triangle there is a 
smaller one (similar to and co-centered with the former) where an equimolar ternary mixture may exist in 
equilibrium with three symmetric phases each rich in one of three components. When the density becomes 
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greater than p q , the equimolar phase is no longer present and the inner four-phase triangle becomes a three- 
phase coexistence region. The three coexisting phases a, [3, and 7 are related symmetrically: phase a is rich 
in component 1 with equal traces of 2 and 3, phase (3 is rich in 2 with equal traces of 1 and 3, and phase 7 is 
rich in 3 with equal traces of 1 and 2. 

Since the three coexisting phases a, [3, and 7 are absolutely symmetric, the three interfacial tensions are 
all equal, a alS = a 13 "' = a ai = a. Therefore a three-phase coexistence will necessarily involve a three 
phase contact region with the contact angles a = (3 = 7 = 2-7r/3. Furthermore, for the appropriate size and 
aspect ratio loq of the cylindrical simulation cell, the point 03 which corresponds to the AL3P partitioning 
will coincide with the phase volume equi-partition point es, and both 03 and es lie exactly in the center of the 
three-phase triangle. 

As an alternative to AL3P, the three-phase partitioning may be of mixed character as illustrated in Fig.4b 
(for 00 > lo' = (6 — 8/7r)/\/3) and Fig.4 (for lo < lo'). Knowing the contact angles (a = (3 = 7 = 27r/3) 
enables one to determine the conditions under which in the state (=03) the AL3P partitioning (Fig.4a) is 
thermodynamically most favorable: 

SH (!<,<,), (30) 



4 - uj aR 

(only the range 1 < lo < it is considered because of condition (10) for the AL2P partitioning). The explicit 
form of the continuous function s(u) is given in the Appendix. The function s(oj)/(4 — u) is positive for 
1 < lo < lo* ~ 2.465, attaining its maximum at lo ~ 1.635. Thus, if r > and 1 < lo < lo*, the AL3P 
partitioning (shown in Fig.4a) is thermodynamically more favorable than those in Fig.4b and Fig.4c. When 
r < (as the mean-field theory predicts 36 for the 3d3c-PS model) and 1 < lo < lo*, making the AL3P 
partitioning thermodynamically most favorable requires that R be sufficiently large, 

4 - to \t\ 



R> R* 



s(oo) a 
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Evaluating R* accurately prior to simulations is complicated by the absence of data on r, even though 
data on a were readily available. Besides, taking large i?'s (for a given density) leads to a significant increase 
in the CPU time for MC simulations. If this is not an acceptable option, one can use a modified version of the 
above method. For the case of a mixed three-phase partitioning, the modifications lead to the equations 

j>(<*) = offMo!)] 1 -* lp(02)] { e- 4HT " T , (31) 

(2 + 7T)UJ 

J_, b(Q2)] g T AuR 1_ 

4R n [p(oi)]«-V(o 3 ) A;T R R 

where a, x, 6, c, d, f are unknown and the explicit form of the function (j)(ui) is given in the Appendix. Gen- 
erating p(mi,m,2), the probability distribution as a function of the order parameters, via appropriate MC 
simulations for a series of i?'s (or L's) at constant go = L/R, and plotting the LHS of eq.(33) vs hiR/R (or 
InL/L), one finds the ratio r/kT as the ordinate intercept of this plot. 

4.2 Monte-Carlo procedure and results 

To obtain the histogram of a system showing a demixing phase transition, most convenient are SGCE MC 
simulations 57-59 , where it is necessary to explore all configurations and compositions of the system for a 
constant total number of molecules. The configurational space is sampled as in a standard canonical simula- 
tion. Varying the composition is carried out by the so-called "identity change" MC moves, and can be done 
in two different ways: (1) choose a molecule at random, change its identity, and accept or reject the change; 
(2) choose a component at random, choose a molecule of that component at random, change its identity, and 
accept or reject the change. These two methods of composition sampling correspond to two different ways 
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to express the SGCE partition function. 57 As a consequence, the acceptance criteria for these methods are 
different. As a slight modification of this algorithm, an identity change attempt may be combined with the 
configurational displacement of a molecule. 60,61 The choice between a pure spatial MC displacement and a 
combined "spatial displacement+identity change" MC move is uniformly random. 

The most obvious advantage of SGCE MC simulations over the Gibbs ensemble simulations 63-65 is 
that less computational effort is required by the former (single phase simulations of N particles) compared 
to the latter (two-phase simulation of N + N molecules). Besides, the Gibbs ensemble method relies on 
particle transfers from one fluid phase (=simulation box) into another, which gives rise to technical limitations 
when studying relatively dense phases, because the probability of particle insertion is then very low and the 
method becomes inefficient. Similar and other problems arise in a grand canonical ensemble simulations of 
dense fluids, although recently some special techniques have been developed to overcome or alleviate those 
problems 23,30 ' 31 ' 51 ' 55 

As a test of the simulation algorithm, the SGCE MC simulations were carried out for the primitive version 
of the two-dimensional Widom-Rowlinson model 56 in a square cell with standard PBC and the critical density 
was determined by using the method proposed by Binder and co-workers. 19 ' 52 ' 53 ' 62 They showed that for 
demixing phase transitions in two-component mixtures, the ratio <s 2 >/<|s|> 2 (the reduced second 
moment of s = 2m — 1, with m being the order parameter) at the critical density p c is independent of the 
system size. Hence, p c in such systems can be determined as the density at which the curves "< s 2 > / < 
|s| > 2 vs density" for different system sizes cross. The simulations were carried out for four different total 
numbers of molecules in the system: N = 256, 726, 1024, and 4098. The results are shown in Figure 
5, where < s 2 > / < \s\ > 2 is plotted as a function of the dimensionless density prj 2 for four iV's. The 
critical density was determined as the average of six densities at which six pairs of simulated series cross. 
The standard deviation of this scatter dominates the standard deviation of < s 2 > / < \s\ > 2 for any given 
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density in a single simulation run. Thus for a two-dimensional Widom-Rowlinson model the dimensionless 
critical density was found to be p c r) 2 ~ 1.58 ± 0.01. This is very close to the result reported previously by 
Johnson et al., 66 who used the invaded cluster (IC) approach in MC simulations: p c r] 2 ~ 1.566 ± 0.003. To 
our knowledge, there have been no other simulations of this model so far. 

The MC simulations of the three-component three-dimensional model were carried out in a cylindrical 
cell with PBC along the z-axis and RPBC at the lateral surface of the cell. Two densities, prj 3 = 0.83 and 
prj 3 = 0.85, were considered, both slightly higher than the quadruple point density of this model reported 
previously 66 as p q r] 3 ~ 0.796. A series of histograms were obtained for each density by carrying out sim- 
ulations for different cylinder sizes L (and R, with uj = L/R = const = 1.5). Equation (25) was then 
applied to each series and the surface tension a for an infinite system was found by extrapolating the series to 
In R/R 2 -► (see Figure 6): arf/kT = 0.0037 ± 0.0005 for prj 3 = 0.83, and arf/kT = 0.0063 ± 0.0008 
for prj 3 = 0.85. As expected, the surface tension increases with increasing density as the three phases be- 
come less and less miscible. Furthermore, according to the mean field approximation, 36 the line tension in 
this model above the quadruple point density is expected to be negative. Hence, prior to obtaining the final 
result, there is some uncertainty whether condition (30) is fulfilled for a given R (with u> = 1.5). The line 
tension was thus determined by using both eq.(21) and eq.(33) and the results were substituted into eq.(30) 
for verification. In neither case condition (30) held whence one can conclude that the dominant three-phase 
partitioning is the one shown in Fig.4b, and eq.(33) is the appropriate one for determining the line tension. 
The results for both densities (as shown in Figure 7) are: rn/kT = —0.093 ± 0.006 for prj 3 = 0.83 and 
TT]/kT = —0.103 ± 0.008 for prj 3 = 0.85. As expected, 36 the line tension is negative. Relatively large 
standard deviations (about 10% of average values) of results for a/kT and r/kT are due to the uncertainty 
in finding the probabilities p(oi), p(o2), and ^(03) from a simulated histogram. To increase the accuracy of 
simulated a/kT and r/kT, simulations longer than 2 x 10 7 MC sweeps (a typical length of one simulation 
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run in the present work) are required. 

5 Conclusions 

A system where the three-phase equilibrium involves a three-phase contact region has been studied. The main 
thermodynamic characteristic of a three-phase contact region is the line tension, defined as the linear excess 
grand canonical potential per unit length of the contact line. In such a system, the probability distribution of 
states of the system with respect to two independent order parameters depends on both the surface and line 
tensions. This constitutes the foundation of a method for determining the line tension by MC simulations 
(either in a semi-grand canonical, or in a grand canonical or in a Gibbs ensemble) proposed in the present 
work. The method is the combination of the analysis of the probability distribution function and finite- 
size scaling, and is the natural development of Binder's approach 18 to studying interfacial phenomena and 
phase equilibria in two-phase systems. MC simulations allow one to obtain the probability distribution of 
states of the system as a normalized histogram. However, the histogram provided by MC simulations and 
results extracted therefrom depend on the size of the simulated system. Carrying out simulations for a series 
of system sizes, and extrapolating the results obtained from the corresponding series of histograms to the 
infinite size, one can determine the line tension of the three-phase contact region and the interfacial tensions 
of all three interfaces in an infinite system (and hence the contact angles). The proposed method is illustrated 
by its application to the three-dimensional ternary fluid mixture, in which molecular pairs of like species do 
not interact whereas those of unlike species interact as hard spheres. The three-phase equilibrium in this 
model above its quadruple-point density involves a three-phase contact region. As expected, the interfacial 
tension increases with increasing density. The behavior of the line tension is also compatible with what 
one would expect based on the mean-field approximation, 36 although the predictions of the latter must be 
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regarded with caution because of the large discrepancy between the mean-field quadruple point density ' 
(p^V ~ 0.601) and the simulated one 66 (p^V ~ 0.796). 

The present work suggests that the most accurate data for the line tension can be obtained by simulations 
in a cylindrical cell under conditions where there is only one three-phase contact line in the cell. For a positive 
line tension this condition can be easily fulfilled by an appropriate choice of the aspect ratio of the cylinder. 
For a negative line tension this condition may require the cylinder to have not only a suitable aspect ratio but 
also a large enough radius which may lead to a drastic increase of simulation time. In this case it may be 
a better choice to conduct simulations under conditions where there are two three-phase contact lines in the 
simulation cell. However, in this case the accuracy of results for the line tension can suffer if the contact lines 
are so close that the inhomogeneities in two three-phase contact regions overlap. In order to avoid this effect 
it is necessary to take longer cylinders which again leads to an increase in simulation time. 

It follows from the foregoing that it is possible to determine the line tension by simulations in a cubic 
cell with PBC in two directions, say x and y, and anti-PBC in the third, say z (the extension of our method 
to this case will be given in a sequel 67 to the present work). Moreover, combining such simulations with the 
canonical ensemble MC simulations, carried out under the same boundary conditions and appropriate choice 
of the order parameters, makes it possible to verify the validity of the new linear adsorption equation 39 
dr = AjA/ij + c a da + cpdfi + c^dj, where dr is the change in the line tension due to the changes 
d\ii (i = 1, 2, 3) of the field variables Hi and Aj is the linear adsorption of component i. The last three 
terms (c a da, etc) on the RHS of the above equation make a difference between the old 2 ' 19 and modified 39 
versions of the linear adsorption equation. At present, there are no recipes to analytically evaluate the co- 
efficients c a , cp, c 7 . However, the proposed method enables one to find dr, while the subsequent canonical 
MC simulations can provide J2i Ajd/Xj. The inequality dr + J2i A«A/ij / could be regarded as an indi- 
rect "simulational" confirmation of the validity of the modified linear adsorption equation (this will also be 
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discussed in the sequel. ). Recently, this inequality was shown to hold in the framework of the square- 
gradient approximation of density-functional theory with a special choice for the local part of the free-energy 
density. Note that in the model system, simulated in the present work, da = d/3 = dj = for any changes 
of field variables. Thus, in this particular model dr = — J2i AjA^j, so another model is needed to simulate 
the validity of the new linear adsorption equation. 

Finally, it should be noted that the proposed method can be easily modified to study wetting phenomena 
and, in particular, to determine the line tension attributed to the region of contact of two equilibrium fluid 
phases with an "inert" rigid substrate. As in the case of two-phase equilibrium discussed in section 2, the 
probability distribution of states of such a system is a function of only one independent order parameter, but 
the two maxima of this function (histogram) will now depend on the corresponding fluid-substrate interfacial 
tensions and hence may not be equal, unlike the situation described in section 2 and Figure 1 . Furthermore, 
the minimum of this function will now depend on the line tension associated with the "fluid-fluid-substrate" 
contact region and may not be located strictly in the middle between the two maxima. A statistical-mechanical 
basis for determining the line tension in such a system by MC simulations will be given in a sequel 67 to the 
present work. 
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Appendix. Explicit form of auxiliary functions s(u) and 4>{u) 



The function s(u) in inequality (30) is denned as 



1 , n 6 4 , , 12 

s(w) = --(3 + 7r)w+-^=7r- — - -^=arccos 5 (w) 

3 



+ ^gH\/«-g 2 n - 2 ^ g ^ [7r-2arccosgH] 2 (1< u < u'), (Al) 

* M = 4 (5 + 27r) + ^ " 37! " i^b < ^ < (A2) 



where is the solution of the equation 



jjvrsH + 1(1 - ^(w))Q(^(w)) = ^-ttlo + 1 (A3) 



with the function Q(x) of some variable x being defined as 



. 3smarccosx— 3xarccosx — 3sin 3 arccosx 

Q x = . (A4) 

1 — x 



The function 0(w) in eq.(32) is defined as 



4>{u) = (2 + tt)u + -^=ir - ~ ^V3arccosg(uj) 

3 



o / O 

+ ^5M\/ a -5 2 M- ^^^y[vr-2arccos 5 (^)] 2 (1< a; < u/), (A5) 

2 13 8 16 

^, = -(2 + .) + ^-^-^ (*/<„<„), (A6> 

with g{uj) defined by eq.(A3). 
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Captions to Figures 1-7 

of the manuscript "HISTOGRAM ANALYSIS AS A METHOD FOR DETERMINING THE LINE TENSION BY 

Monte-Carlo simulataions" by Yuri Djikaev. 

Figure 1: Typical "histogram", probability distribution p of states of a system with respect to a single order 
parameter m, under conditions of two-phase coexistence. The values of the order parameter in two equilib- 
rium phases are m a and mr 

Figure 2: Typical probability distribution p of states of a system with respect to two order parameters mi and 
771,2 under conditions of three-phase coexistence, obtained for a three-dimensional ternary fluid mixture, in 
which molecular pairs of like species do not interact whereas those of unlike species interact as hard spheres, 
at the total density prj 3 = 0.85. a) Three-dimensional plot; b) Projection of the histogram onto the mole- 
fraction plane (for details see the text). 

Figure 3: Three two-phase (half-and-half) partitionings compatible with the periodic boundary conditions 
along the z axis: a) AL2P partitioning; b) droplet-like partitioning; c) transversal partitioning. 

Figure 4: a) The axial longitudinal three-phase (AL3P) partitioning in the cylindrical simulation cell with the 
contact angles a, [3, and 7; b) and c) - the alternative to AL3P three-phase partitionings in the cilinder for 

lo > uj' and uj < uj' respectively 

Figure 5: Determination of the critical density for the primitive version of the two-dimensional Widom- 
Rowlinson model by SGCE MC simulations. The reduced second moment (s 2 )/(\s\) 2 (s = 2m — 1 with 
m the mole fraction of one of two components) is plotted vs dimensionless density prj 2 for four different 
numbers of molecules in the simulated system (square with PBC): larger circles are for N = 256, smaller 
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points for N = 726, smaller circles for N = 1024, and larger points for N = 4098. The thin curves are 
plotted only for guiding the eye. The dashed vertical line indicates the critical density from ref.66. 

Figure 6: Determination of the surface tension by the histogram analysis combined with the finite size scal- 
ing. The LHS of eq.(25) is plotted vs InL/L 2 for two different densities, pif = 0.83 (upper series) and 
prf 3 = 0.85 (lower series). The straight lines show the linear least-squares fits to each series of simulated 
points. The extrapolated values for the infinite system size are: arj 2 jkT = 0.0037 ± 0.0005 for prf = 0.83, 
and aif/kT = 0.0063 ± 0.0008 for P 7] 3 = 0.85. 

Figure 7: Determination of the surface tension by the histogram analysis combined with the finite size scal- 
ing. The LHS of eq.(33) is plotted vs In L/L for two different densities, prj 3 = 0.83 (upper series) and 
prf = 0.85 (lower series). The straight lines show the linear least-squares fits to each series of simulated 
points. The extrapolated values for the infinite system size are: T-q/kT = —0.093 ± 0.006 for rhor) 3 = 0.83, 
and Trj/kT = -0.103 ± 0.008 for rhor] 3 = 0.83. 
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